#  Table b2
rm(list=ls())

(WD <- getwd())
if (!is.null(WD)) setwd(WD)

indata0      = paste0(WD,"/_individual_data/_spfrawdata/","", collapse = NULL)
indata1      = paste0(WD,"/_other_data/_public_signals/_hard_variables/","", collapse = NULL)
indata2      = paste0(WD,"/_other_data/_public_signals/_survey_variables/","", collapse = NULL)

library(haven)
library(AER)
library(sandwich)
library(lmtest)
library(pracma)
library(stargazer)
library(plm)
library(pracma)
library(DataCombine)
library(jtools)
library(plyr)
library(ggplot2)
library(tidyverse)
library(dotwhisker)

lagpad <- function(x, k) {
  res <- c(rep(NA, k), x)[1:length(x)]
  return(res)
}

n_hard    = 9
n_survey  = 4

load(paste(indata0,"us_spf_pgdp_ind.Rda",sep=""))
data$ID         = data$id
data            = data[data$qdate>=1970.00,]
data            = data[data$qdate<=2020.00,]
data            = ddply(data,.(qdate),transform, cons = mean(f2,na.rm = TRUE))
data_fcast      = subset(data, select = c(qdate,fc_err, ID, p_realz))

data_tmp        = aggregate(. ~ qdate, data_fcast, mean)
data_tmp$lag    = lagpad(data_tmp$p_realz,5)
data_tmp        = subset(data_tmp, select = c(qdate, lag))

load(paste(indata1,"public_signals_hard.Rda",sep=""))
data_hard       = data
data_hard$qdate = data_hard$qdate+5/4 
data_hard       = merge(data_hard, data_fcast, by = "qdate", all = TRUE)
data_hard       = merge(data_hard, data_tmp, by= "qdate", all = TRUE)

load(paste(indata2,"public_signals_survey.Rda",sep=""))
data_survey       = data
data_survey$qdate = data_survey$qdate+5/4
data_survey       = merge(data_survey, data_fcast, by = "qdate", all = TRUE)

data_all       = merge(data_hard, data_survey, all = TRUE)

data_all$z_dlneer = -scale(data_all$dlneer)
data_all$z_dlimp  = +scale(data_all$dlimp)
data_all$z_dloil  = +scale(data_all$dloil)
data_all$z_cf_exp_infl = +scale(data_all$cf_exp_infl)
data_all$z_dlstocks  = -scale(data_all$dlstocks)
data_all$z_term  = -scale(data_all$term)
data_all$z_tips  = +scale(data_all$tips)
data_all$z_urate  = -scale(data_all$urate)
data_all$z_lag  = -scale(data_all$lag)
data_all$z_spf  = +data_all$spf
data_all$z_liv  = +data_all$liv
data_all$z_mich = +data_all$mich
data_all$z_sce  = -data_all$sce

plm_all     = plm(fc_err ~ z_dlneer + z_dlimp + z_dloil + z_cf_exp_infl+ z_dlstocks + z_term + z_urate + z_lag + z_spf, data=data_all,index=c("ID", "qdate"), model="within")
beta        = as.numeric(coefficients(plm_all))
std_error   = as.numeric(sqrt(diag(vcovDC(plm_all, type = "HC1"))))

multi_est  =  data.frame(beta, std_error)

term = c('NEER','IMPP','OIL','FIN', 'STOX','TERM','UNMP', 'LAG','SPF')
multi_est$term = term
names(multi_est)[1] <- "estimate"
names(multi_est)[2] <- "std.error"
multi_est = multi_est[c(9,8,1,2,3,7,4,5,6),]

write.table(format(multi_est,digits=3), file = "_tables/table_b2_top.txt", sep = "\t")

plm_all     = plm(fc_err ~ z_spf + z_sce +z_mich + z_liv, data=data_all,index=c("ID", "qdate"), model="within")
beta        = as.numeric(coefficients(plm_all))
std_error   = as.numeric(sqrt(diag(vcovDC(plm_all, type = "HC1"))))

multi_est  =  data.frame(beta, std_error)

term = c('SPF','SCE','MICH','LIV')
multi_est$term = term
names(multi_est)[1] <- "estimate"
names(multi_est)[2] <- "std.error"
multi_est = multi_est[c(1,3,2,4),]

write.table(format(multi_est,digits=3), file = "_tables/table_b2_bottom.txt", sep = "\t")

